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The dynamics of samples in the continuous-imaginary-time quantum world-line Monte Carlo 
simulations with extended ensembles are investigated. In the case of a conventional flat ensemble 
on the one-dimensional quantum 5 = 1 bi-quadratic model, the asymmetric behavior of Monte 
Carlo samples appears in the diffusion process in the space of the number of vertices. We 
prove that a local diffusivity is asymptotically proportional to the number of vertices, and we 
demonstrate the asymmetric behavior in the flat ensemble case. On the basis of the asymptotic 
form, we propose the weight of an optimal ensemble as 1/ ^/n, where n denotes the number of 
vertices in a sample. It is shown that the asymmetric behavior completely vanishes in the case 
of the proposed ensemble on the one-dimensional quantum 5* = 1 bi-quadratic model. 
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Properties of various quantum exotic states and phase 
transitions between them have been extensively investi- 
gated. For example, quantum paramagnetic and valence- 
bond-soUd states related to a possible mechanism to sup- 
port a novel superconductivity in cuprates,^-* and the ex- 
istence of non-Landau-Ginzburg- Wilson type phase tran- 
sitions between different broken symmetries. In or- 
der to numerically investigate such quantum states and 
phenomena, unbiased quantum world-line Monte Carlo 
(QMC) methods based on a Markov chain are powerful 
tools, because they can be applied to large scale sys- 
tems at low temperatures and are not limited to the 
one-dimcnsional case, if no negative-sign problem exists. 
However, even for no negative-sign cases, it is sometimes 
difficult to get accurate data from the conventional QMC 
simulations, because Monte Carlo samples are trapped 
near a metastable state in configuration space. In general, 
the quality of QMC simulations deteriorates at low tem- 
peratures, because metastable states are related to de- 
generate ground states. In fact, the autocorrelation time 
of samples in QMC simulations becomes exponentially 
large to overcome barriers between metastable states. In 
order to eliminate the rapid increase in autocorrelation 
times in Monte Carlo simulations, two approaches were 
proposed over the last two decades. One is to change the 
local update method of samples in a Markov process to 
a global one that connects metastable states directly. In 
fact, the loop algorithm^^ has been successful in studies 
of quantum magnetic phases due to the global update 
whose shape is like a loop and which corresponds to a 
magnetic correlated domain. However, in some cases us- 
ing the loop algorithm, we encountered rapid increases 
in autocorrelation time: For example, valence-bond-solid 
states that break spatial symmetries on the two and 
quasi-one dimensional lattice.'*'^-' For such cases, the sec- 
ond approach is probably effective, in which a canonical 
ensemble is replaced by an artificial extended ensemble 
such that Monte Carlo samples would not be trapped 
near metastable states. For classical systems, the ex- 
tended weight is adjusted such that the appearance ratio 



of energy E samples would be fiat. They are not trapped 
near a metastable state because Monte Carlo samples 
diffuse in a wide energy range. Therefore, these meth- 
ods have been extensively used for studies of spin glasses 
and frustrated classical models. However, for quantum 
models, they have been tested only in a few cases. ^'^^ In 
order to get conclusive numerical results for quantum 
strong correlated phenomena as mentioned above, we 
need to understand the property of QMC methods with 
extended ensembles, and it is important to improve their 
efficiencies. In this letter, we will concentrate the dif- 
fusive behavior of samples in the continuous-imaginary- 
time QMC simulations with an extended ensemble. We 
will report the asymmetric behavior in the diffusion pro- 
cess for the one-dimensional quantum 5=1 bi-quadratic 
(BQ) model case. On the basis of our proven asymptotic 
form of the local diffusivity of samples, we will propose an 
optimal ensemble. The performance of the proposed en- 
semble will be shown for the case of the one-dimensional 
BQ model. 

Although the first formulation of QMC methods 
with extended ensembles (EEQMC) was based on high- 
temperature series expansion,^' it can be also done 
on the path-integral representation with a continuous- 
imaginary-time limit (see §2.15 in rcf. 8). In order to 
provide a brief description of the EEQMC algorithms 
on the path-integral representation, we first consider the 
definition of the exponential operator: 



exp 



lim 



n 



1 - 



M 



(1) 



Inserting the identity operator 1 = X]c(l"^)('^l with a 
complete orthonormal basis {|a)} between two adjacent 
factors in the right-hand side of eq. (1), we obtain the 
discrete-imaginary-time path-integral representation of a 
partition function as 

M K 
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where Tib is the 6-th interaction Hamihonian, (3 is an 
inverse temperature, A = 1/M, SK+i{i) = Si{t + 1), 
and 5*1 (M + 1) = 5*1(1). Next, we introduce new auxihary 
variables Gb{t) cahed graph variables as (1 — /37ibA) = 
EG6(t)=o.i(-/3^f'^)'^''^*^- Then, eq. (2) is rewritten as 

S G n 

Wo{s,G) = \{{SuM-nu^f- \Su), (4) 

u 

^N = E E ^o(5,G), (5) 

S (G|ri(G)=n) 

where u, w', 5, and G denote (&, t), (6 + {5'„}, and 
{Gu}, respectively, and n{G) = X^m^m- I'^ the follow- 
ing, a graph variable G„ that takes a value 1 is called 
a vertex and the number of vertices n[G) is called a 
vertex number. This representation is the mathematical 
background required to describe the remarkable QMC 
algorithms that have been developed during the last 
two decades.®^ In particular, we can take a limit of a 
continuous-imaginary time. A 0, in the level of QMC 
algorithms (see §2.5 in ref. 8 for details). In the remain- 
der of this paper, we consider the EEQMC algorithm on 
the continuous-imaginary time, because it has no system- 
atic error from the discretization of an imaginary time. 

The vertex number n(G) in a canonical ensemble sta- 
tistically corresponds to an inverse temperature (i be- 
cause {n{G))p = P(—H)i3, where {■)/3 is a canonical en- 
semble average. Therefore, if we adjust the weight of 
a configuration (5*, G) so that the frequency of obtain- 
ing the vertex number n would be independent of n, 
i.e., flat,^' we can sample various configurations in a 
wide inverse temperature range. From eq. (3), n{n) is 
regarded as the density of states with a fixed vertex num- 
ber n. Hence, if the factor /S"^'^' in eq. (3) is replaced by 
l/J7(n(G)), it can be done. In general, in order that the 
appearance ratio of configurations with a vertex number 
n is Py (n) , the extended ensemble weight of a configura- 
tion (5,G) has to be Pyin{G))Wo{S,G)/n{n{G)). How- 
ever, we need to guess ^l{n) from the QMC samples them- 
selves, because f2(n) is not known a priori. Fortunately, 
some sophisticated methods have been proposed.^' In 
our simulations, we have used the broad-histogram rela- 
tion for the vertex number^^^ as 

n{n + l) ^ (diag(-H))„ 
Q{n) n + 1 - {nK)n+i ' 

where {Q)n denotes the micro-canonical ensemble aver- 
age of an operator Q at a fixed vertex number n, diag((5) 
refers to the diagonal part of an operator Q, and hk 
is the number of kinks at which a state changes, i.e., 
{Su'\Su) = 0. Because the right-hand side in eq. (6) can 
be directly estimated in the EEQMC simulations, n{n) 
can be calculated from this recursion formula. We should 
note that this estimation method is independent of the 
dynamics of the EEQMC samples; it is not based on the 
histogram of the appearance of a vertex number in the 
EEQMC simulations. 

However, after ^l{n) is sufficiently estimated, the dy- 



namics of the EEQMC samples do not seem to be a reg- 
ular random walk in the vertex number space. In par- 
ticular, in order to investigate this behavior, we focus 
on the first-passage time (FPT) of EEQMC samples re- 
garded as random walkers in a vertex number space. The 
FPT is defined as the time at which a random walker 
first reaches a threshold value. Because the movement of 
an EEQMC sample is usually restricted to the interval 
[Na, Nf,] in the vertex number space, two types of FPTs 
are defined as a QMC sample moves to Nf, from Na, and 
vice versa. In the following, the former is called forward 
and the latter is called backward. 
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Fig. 1. First-passage times (FPT) of EEQMC samples for the 
one-dimensional quantum 5 = 1 bi-quadratic model. Filled tri- 
angle and square symbols are forward and backward FPTs for a 
flat ensemble, respectively. Triangle and square symbols are for- 
ward and backward FPTs for the ensemble defined by cq. (15), 
respectively. In the inset, the ratios of backward and forward 
FPTs are shown. 

Figure 1 shows forward and backward FPTs, Tfp(0 
Nm) and rFp(A^rn 0), for the one-dimensional BQ 
model with chain lengths L = 8,16,32,64,128, and 
256. The Hamiltonian of the BQ model on the one- 
dimensional lattice is H = ' Si+i)^, where Si 
denotes an S* = 1 spin operator at the i-th site. Because 
this model is integrable, the ground state is well known as 
the dimerized state in which singlet pairs align and com- 
pletely cover the one-dimensional lattice. This dimerized 
state is twofold degenerate and spontaneously breaks a 
translational symmetry. And the autocorrelation time in 
QMC simulations is the round-trip time between such 
degenerate states. As mentioned above, even if we use 
a loop algorithm for the BQ model, ^-^^ it grows rapidly 
at low temperatures. Therefore, this model is suited to 
the test of the EEQMC algorithms. In all cases, is 
proportional to the chain length L as Nm = 64i. It cor- 
responds to a constant low temperature enough to cal- 
culate a ground state, because the one-dimensional BQ 
model has a gap. The total number of EEQMC sam- 
ples for a chain length is sufficiently large for estimating 
FPTs: For example, 32 independent runs with 6 x 10^ 
Monte Carlo sweeps (MCSs) for L = 256. We should 
note that MCS is adopted as a unit of time. One MCS for 
a configuration (S*, G) in the continuous-imaginary-time 
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EEQMC algorithm consists of three steps: (i) deciding 
a new vertex number n' under a given S'-configuration, 
(ii) assigning new graphs G' with the vertex number n' 
to a given S'-configuration, and (iii) choosing a new S'- 
configuration under the given graphs G' . It is possible 
to measure observables only at the time that these three 
steps are completed. In Fig. 1, forward and backward 
FPTs for a flat ensemble increase almost linearly, but 
the backward FPT is always larger than the forward one: 
For example, the ratio of two FPTs is 7.0(5) for L ~ 256 
(see the inset of Fig. 1). Thus, the EEQMC samples for 
a flat ensemble move quickly from high temperatures to 
low ones, but slowly in the reverse direction. In order to 
improve the efficiency of the EEQMC algorithms, it is 
necessary to correct this asymmetric behavior. 

When we make a new G' configuration under a fixed 
S'-configuration, the vertex (G„ = 1) at a kink can not be 
removed, because the local weights in Wq{S, G) at a kink 
becomes zero: {Su'\Su) ~ 0. Therefore, if the number of 
kinks is nx^ the probability p{n'\nK) to choose the next 
vertex number n' in step (i) is proportional to the sum 
of weights of the configurations that have unchanged nx 
vertices at kinks and new inserted [n' — nx) ones into a 
given S'-configuration: 

-1 {n' — n^ 

Pv{n') E(.|G„=o)(^«'l(-W«A)|S„) 



p{n'\nx) (X 



(7) 

For finite-size systems, if the vertex number is sufficiently 
large, the right-hand side in eq. (6) is converged. Using 
limiting values as wq = lim„^oo (diag(— 7Y))^ and rx ~ 
lim„_>oo {nx)n l^-i asymptotic form of f2(n) is as 



17 (71) 



VL(n- 1) 



Wo 



1 - rx 



1 

n! 



Wo 



\^rx 



(8) 



Substituting eq. (8) into eq. (7), we find that the main 
factor of 'p(ri!\nx) is the negative binomial distribu- 
tion. Using the limit theorem for the negative bino- 
mial distribution, we obtain the asymptotic form of 
■p{Ti!\nx) oc P„(n') exp[— (n' — mo)^/2(mo/rif )], where 
mo = {r~^ — l)nx- Here, we assume that Pv{n') is a 
slowly varying function. And if we assume that the prob- 
ability p{nx\n) that the number of kinks in a configura- 
tion with a vertex number n is nx is equivalent to the 
probability p(rl\nx)^ the probability to choose the next 
vertex number n' from configurations with a vertex num- 
ber n is 



p{n^n')= I dnx p{n'\nx)p{nx\n), 





'2{l-rx)' 











(9) 



(10) 



where NQ{m, cr^) denotes the Gaussian distribution with 
mean m and variance cr^. Thus, the EEQMC samples 
almost seem to be random walkers in the vertex number 
space, but the local diffusivity D{ti) increases linearly as 



D{n) 



(11) 



ure 2 shows the local diffusivity D{n) for a flat ensem- 
ble (Pv{n) = 1) in the EEQMC simulations of the one- 
dimensional BQ model. The chain length L is 128 and the 
total number of MCSs is 1.05 x lO'^. The local diffusivity 
in Fig. 2 is approximately linear in the region above the 
vertex number 100 (see also the left-top inset of Fig. 2). 
The solid line in Fig. 2 is a linear function predicted in 
eq. (11) with rx — 0.6108(1), which is evaluated from 
the QMC simulations. The predicted line is consistent 
with the local diffusivity in the region above the ver- 
tex number 1000. And the discrepancy between them is 
never more than 6% at all vertex numbers below 1000 
but zero. Next, we checked the dependence of the local 
diffusivity on extended ensembles. In the right-bottom 
inset of Fig. 2, the ratios between local diffusivitics in 
two different ensembles are shown. Because these values 
are approximately equal to one, the local diffusivity is al- 
most unaffected by the choice of extended ensembles. For 
other system-size cases, the same results were obtained. 
Therefore, the local diffusivity of the EEQMC samples 
is described well by eq. (11). 
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Fig. 2. Local diffusivity in the EEQMC simulations of the one- 
dimensional quantum 5 = 1 bi-quadratic model. The chain 
length L is 128. The extended ensemble is flat: Pvin) = 1. The 
solid line shows the predicted linear increase in eq. (11). The 
value of tk is 0.6108(1), which is evaluated from the QMC sim- 
ulations. In the left-top inset, the local diffusivity in the small 
vertex number region is shown. In the right-bottom inset, the ra- 
tios between local diffusivitics in a flat ensemble and the optimal 
one defined by eq. (15) are shown. 



From eq. (10), the behavior of the samples in 
the EEQMC simulations may be described well by 
a Fokker-Planck equation (FPE) on a vertex number 
range [Na,Ni]}^"> Using the theory of ffist -passage pro- 
cesses,^^' we can explicitly obtain the first-passage times 
for a one-dimensional FPE as 



Tfp{Na ^ n) 



Tfp{Nb ^ n) 



dxPv (x) 



dx' 



dxPy (x) 



D{x')P^{x'y 
dx' 



(12) 



(13) 



where n is the number of vertices in a configuration. Fig- 



D{x')P^{x') 

If we assume the linear increase of the local diffusivity 
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D{ti) as in cq. (11), FPTs for a flat ensemble {Pv{n) ~ 1) 
are as 

iVn, TFLA^(iV^ ^ 0) 



^FP 



(0 ^ N^) 



' K 



1 



logiV^ 



(14) 

These results are qualitatively consistent with the behav- 
ior of FPTs in Fig. 1. The forward FPT increases linearly 
and the backward one is always larger than the forward 
one and their ratio varies slowly. Thus, the qualitative 
behavior of the EEQMC samples is described well by 
the FPE with eq. (11). 

In order to correct the asymmetry between forward 
and backward FPTs for a flat ensemble, an extended en- 
semble Pv{n) should be adjusted so that the right-hand 
side in eq. (12) would be equivalent to that in eq. (13). 
While it typically cannot be uniquely determined, a spe- 
cial solution exists: Pv{n) = \/ ^ D{n). This special en- 
semble was first derived from a maximization of random 
walker flows. Using this ensemble, we find that the 
forward FPT is not only equivalent to the backward one, 
but also the total FPT, Tfp(0 ^ N^) + Tfp{N^ 0), 
is minimized. Therefore, from eq. (11), the optimal 
ensemble for the EEQMC methods is 

Pnn) = -ji= cx ^. (15) 

Substituting eqs. (11) and (15) into eqs. (12) and (13), we 
find that the forward and backward FPTs for the optimal 
ensemble are equivalent and become twice as large as the 
forward one for a flat ensemble: 



2rFLAT(0_>7V^). 



Tp"p^'(7v.^o) = rF"p^'(o 

(16) 

The forward and backward FPTs in the EEQMC simu- 
lations with the optimal ensemble are shown in Fig. 1. 
Here, although eq. (11) is the asymptotic form, we use 
eq. (15) for all vertex numbers except zero, and the value 
at the vertex number zero is defined by that at the ver- 
tex number one: Pv{0) = P„(l). In Fig. 1, the forward 
FPT for our proposed ensemble is consistent with the 
backward one (see also the inset of Fig. 1). And, as in 
eq. (16), they are approximately twice as large as the 
forward FPT for a flat ensemble. Thus, the asymmetry 
is completely corrected by the ensemble in eq. (15). 

As mentioned above, we adopt MCS as the unit of 
time. But the actual computational times for one MCS 
are not constant. In fact, in many cases the number of 
local steps in one MCS is proportional to the vertex num- 
ber in a sample. Therefore, the local step can be adopted 
as the unit of time. In this case, from eq. (11), we find 
that the local diffusivity in units of local steps is indepen- 
dent of the vertex number. Hence, the optimal ensemble 
in units of local steps is flat. In other words, the optimal 
ensemble for computational times is as P^'^^^(n) = 1/n. 
The total local steps of FPTs is | times that for P°^"^. 
But the number of samples in the region of large vertex 
numbers decreases more than that for Py^^- In order to 
calculate a canonical ensemble average of an observable 
at an inverse temperature (3, it is necessary to calculate 



the reweighted summation of micro-canonical ensemble 
averages. Vertex numbers that mainly contribute to it 
are in the region of which the width is proportional to 
yjn{(3), where n{j3) is a center of the region. Therefore, 
in the case of P^^"^ , the number of samples that con- 
tribute to a canonical ensemble average is constant at all 
inverse temperatures. But that for the ensemble Py''^'^ 
decreases at low temperatures. 

In summary, we considered the diffusion of samples in 
the continuous-imaginary-time EEQMC simulations. In 
particular, the asymmetric behavior of FPTs of EEQMC 
samples was reported in detail. We proved that the local 
diffusivity of the EEQMC samples is asymptotically pro- 
portional to the vertex number. And it was shown that 
the asymptotic form is consistent with the local diffu- 
sivity in the EEQMC simulations of the one-dimensional 
BQ model in the wide region of the vertex numbers. Us- 
ing this result and the theory of first-passage processes, 
we demonstrated the asymmetric behavior of FPTs for a 
fiat ensemble case and proposed an optimal ensemble for 
the continuous-imaginary-timc EEQMC simulations in 
order to correct the asymmetric behavior. It was shown 
that the asymmetric behavior on the one-dimensional 
BQ model completely vanishes in the case of the pro- 
posed ensemble. 
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